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Abstract We study the nature of motion in a 3D po- 
tential composed of perturbed elliptic oscillators. Our 
technique is to use the results obtained from the 2D po- 
tential in order to find the initial conditions generating 
regular or chaotic orbits in the 3D potential. Both 2D 
and 3D potentials display exact periodic orbits together 
with extended chaotic regions. Numerical experiments 
suggest, that the degree of chaos increases rapidly, as 
the energy of the test particle increases. About 97% of 
the phase plane of the 2D system is covered by chaotic 
orbits for large energies. The regular or chaotic charac- 
ter of the 2D orbits is checked using the 5(c) dynamical 
spectrum, while for the 3D potential we use the S(c) 
spectrum, along with the P(f) spectral method. Com- 
parison with other dynamical indicators shows that the 
S(c) spectrum gives fast and reliable information about 
the character of motion. 

Keywords Galaxies: kinematics and dynamics; 
dynamical indicators; elliptic oscillators 



1 Introduction 

It has been almost five decades, since the pioneer work 
of Hcnon and Heiles [1], who studied the regular and 
chaotic character of motion in a system of two cou- 
pled harmonic oscillators. During these years several 
methods have been presented in order to distinguish 
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order from chaos. Among them, it is worth mention- 
ing, the method based on the Fourier spectra of trajec- 
tories (see [2,18]), the Lyapunov Characteristic Expo- 
nents (L.C.Es)(see [19]) and the KS entropy (see [3]). 

Apart from the above methods, in systems of two de- 
grees of freedom, one can apply the classical method of 
the Poincare phase plane, in order to distinguish regular 
from chaotic motion. On the contrary, in a 3D system 
the above method can not be applied. On this basis, 
it is of particular interest to locate the regions of ini- 
tial conditions in 3D potentials generating regular or 
chaotic motion. 

The aim of the present article is to investigate the 
nature of motion in the 3D potential 

,.,2 

V{x,y,z) 



(x 2 + y 2 + z 2 ) 



u> 
~2 

+ e [x 2 y 2 + y 2 z 2 + x 2 z 2 - x 2 y 2 z 2 ] , (1) 

where u> is the common frequency of oscillation along 
x, y and z axis, while e is the perturbation strength. Po- 
tential (1) represents three coupled harmonic oscillators 
in the 1:1:1 resonance. Potentials of this type are also 
known as perturbed elliptic oscillators (see [5,6]). The 
basic reason for the choice of potential (1) is that per- 
turbed elliptic oscillators appear very often in galactic 
dynamics and atomic-particle physics (see [5] and refer- 
ences therein). A second reason is that it displays exact 
periodic orbits, interesting sticky orbits together with 
large chaotic regions. Therefore, it gives a good chance 
in order to check the reliability of the S(c) spectrum 
(see [7,10]), which is used as an indicator to track the 
sticky orbits and to distinguish between regular and 
chaotic motion. 

The outcomes of this research are mainly based on 
the numerical integration of the equations of motion 
d V(x, y, z) 



x = — - 



dx 
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dV(x,y,z) 

y = s — ' 

dy 

= d V(x,y,z) 
dz 

(2) 

where the dot indicates derivatives with respect to the 
time. 

The Hamiltonian to the potential (1) reads 
H=^{pl+p 2 y +p 2 z )+V(x,y,z) = h , (3) 

where p x ,p y and p z are the momenta per unit mass 
conjugate to x, y and z, while h is the numerical value 
of the Hamiltonian. 

The paper is organized as follows: In Section 2 we 
present the analysis of the structure of the (x,p x ) Poincare 
phase plane of the 2D system for a set of values of en- 
ergy h. In the same Section, we give some examples 
for the reader to understand the usefulness of the S(c) 
spectrum, since it has the ability to describe islandic 
and sticky motion. In Section 3 we study the character 
of orbits in the 3D system. Special interest is given to 
the evolution of the exact periodic orbits. Particularly, 
we are interested to find the value of energy, where or- 
bits starting near the exact periodic orbits change their 
character from regular to chaotic. We close our investi- 
gation with Section 4, where a discussion and the con- 
clusions of this research are presented. 

2 Order and chaos in the 2D system 

In the following we shall study the character of motion 
in the corresponding 2D potential 

V(x,y) = ^(x 2 +y 2 )+exV • (4) 

It is well known, that in this potential the two axes x 
and y and the straight lines x = ±y are exact periodic 
orbits (see [8] ) . The Hamiltonian to potential (4) writes 

H 2 = 1 -(pl+p 2 y )+V(x,y) 

= \(PI+P 2 y+ u 2 x 2 + w V) + ec V 

= h 2 , (5) 

where h 2 is the numerical value of the Hamiltonian. 
It is interesting to note, that potential (4) has not a 
finite energy of escape, so the corresponding zero veloc- 
ity curves (ZVCs) are always closed. In all numerical 
calculations, we adopt the values: to = l,e = 1, while 
the energy is treated as a parameter. 

Figure la shows the (x,p x ), (y = 0,p y > 0) phase 
plane, when h 2 = 1. One can see, that a large part of 



the phase plane is covered by chaotic orbits. The reg- 
ular regions are occupied by invariant curves around 
the points (x ,p x0 ) = (0,0) and (x ,p x o) = (o,±\/Jij. 
The above points give the positions of two exact pe- 
riodic orbits, The first is the y axis while the second 
describes the x — ±y straight line orbits on the (x,p x ) 
phase plane. Furthermore, one observes a large num- 
ber of smaller islands of invariant curves produced by 
secondary resonances and some sticky regions as well. 
Figure lb is similar to Fig. la but when h 2 = 2. Here 
we see that the majority of the phase plane is covered 
by a chaotic sea. The regular regions are confined only 
around the straight line periodic orbits, while the y 
axis has now become unstable. Some smaller islands 
of invariant curves are also present. The most interest- 
ing feature observed in this case, is the sticky regions 
around the chain of islands of invariant curves on the x 
as well as on the p x axis. We shall come to this point 
later in this Section. If we increase the value of energy to 
h 2 = 4, we obtain the phase plane shown in Figure lc. 
In this case, the y axis has returned to stability, while 
the straight line periodic orbits are unstable. Almost 
all the phase plane is chaotic, except of a small region 
around the center and some small islands of invariant 
curves which are products of secondary resonances. A 
sticky region around the chain of the small islands of 
invariant curves near the center is also observed. Fig- 
ure Id is similar to Fig. la but when h 2 = 7.174235. 
This value of the energy is the value of the energy of 
escape of the 3D potential (see next Section). Here al- 
most all the phase plane is covered by chaotic orbits. 
There is a tiny regular region near the center. A careful 
observation also shows some very small regular islands 
of invariant curves embedded in the chaotic sea. 

One can see, that the fraction of the phase plane 
covered by chaotic orbits increases as the energy h 2 in- 
creases. Figure 2a shows a plot of the percentage of 
the phase plane A% covered by chaotic orbits vs h 2 . 
We observe that A% increases rapidly, as h increases. 
Dots indicate values of A% found numerically, while the 
solid line is a fourth degree polynomial fitting curve. 
We would like to make clear, that A% is estimated on 
a completely empirical basis, by measuring the area on 
the (x,p x ) phase plane occupied by chaotic orbits (see 
[11]). In order to have an estimation of the degree of 
chaos from another point of view, we have plotted the 
maximum L.C.E vs h 2 . The results are shown in Fig- 
ure 2b. Note that the L.C.E increases linearly with h 2 . 
Here, we must point out that it is well known that the 
L.C.E has different values in each chaotic component 
(see [20]). Since we have regular regions in all cases and 
a large chaotic area with small sticky regions embedded 
in each chaotic sea, we calculate the average value of the 
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Fig. 1 (a-d): The (x,p x ) phase plane, when ui = l,e = 1, while (a-upper left): h 2 = 1, (b-upper right): h 2 = 2, (c-down left): 
h 2 = 4 and (d-doum right): h 2 = 7.174235. 





Fig. 2 (a-b): (a-left) A plot of A% vs 7i 2 and (b-right) a plot of L.C.E vs h 2 . 



L.C.E in each case by taking fifty orbits with different Interesting information not only for the chaotic or 

initial conditions (xo,p x o), in each chaotic domain. In regular character of orbits but also for islandic and 
all cases, the obtained values of the L.C.Es were dif- sticky motion can be obtained using the dynamical pa- 
ferent in the fourth decimal point, in the same chaotic rameter of the S(c) spectrum. The parameter c is de- 
area, fined as 

Pxi 



Ci = 



Pyi 



(6) 
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Fig. 5 (a-d): (a-upper left): The S(c) spectrum of the regular orbit producing the two sets of five islands of invariant curves 
in the (x,p x ) phase plane, (b-d): The evolution of the 5(c) spectrum of a sticky orbit when (b-upper right): T = 1000, (c-down 
left): T = 4200 and (d-down right): T = 4300. 
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where {xi,p X i,p y i) are the successive values of the (x,p x , 
orbital elements on the Poincare (x,p x ) phase plane. 
The dynamical spectrum of the parameter c is its dis- 
tribution function 

AN{c) 



Py) 



5(c) = 



NAc 



(7) 



where AN(c) are the number of the parameters c in the 
interval (c, c+ Ac), after TV iterations. This dynamical 
spectrum has been proved a very reliable tool in sev- 
eral cases (see [7,10] for details). The nature of a 2D 
orbit can be identified by the shape of the 5(c) spec- 
trum. If the structure of the spectrum is a well formed 
i7-type shape, then the corresponding orbit is regular. 
On the other hand, if the shape of the spectrum is 
complicated, asymmetric, with a lot of of small and 
large abrupt peaks, then this indicates that the orbit is 
chaotic. Moreover, the 5(c) spectrum can help us iden- 
tify resonant orbits or orbits of higher multiplicity, as it 
produces as much /7-type spectra as the total number 
of the islands of invariant curves in the (x,p x ) phase 
plane. One additional advantage of this spectrum, is 
that it can be used in order to calculate the sticky pe- 
riod of a 2D sticky orbit and also to help us follow its 
time evolution towards the chaotic sea. 

At this point, we should emphasize that there are 
some spectra, such as the 5(a) spectrum, (see [23]), 
which have no ability to identify islandic or sticky mo- 
tion. Furthermore, the 5(a) spectrum needs the com- 
putation of two neighboring orbits and a large number 
of iterations (about 5 x 10 4 ) in order to provide reli- 
able results regarding the character of a 2D orbit. Us- 
ing the 5(c) spectrum we have managed to overcome 
these drawbacks. Recently, in [24] we introduced the 
S(g) spectrum which proved to be a very sensitive and 
reliable spectrum. This new spectrum definition allow 
us to study more accurately the islandic motion of reso- 
nant orbits of higher multiplicity, as it produces as much 
spectra as the total number of the islands of invariant 
curves in the (x,p x ) phase plane. The main advantage 
of the S(g) is that it has the ability to detect com- 
plicated resonant orbits especially in 2D Hamiltonian 
systems, no matter how small are the corresponding is- 
lands of invariant curves. Furthermore, for the study of 
orbits in 3D Hamiltonian systems, we constructed and 
used in [25] the S(w) dynamical spectrum. This spec- 
trum is an advanced form of the 5(c) spectrum and it 
was especially designed for the study of 3D orbits, as 
it produces as much spectra as the total number of the 
3D invariant tori in the (x,p x ,z) phase subspace. 

In what follows we present some examples in order 
to visualize the usefulness of the 5(c) spectrum. Fig- 
ure 3a-b shows the 5(c) spectrum and the correspond- 
ing L.C.E vs time of an orbit with initial conditions: 
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Fig. 6 The L.C.E vs time for the sticky orbit of Fig. 5. 




Fig. 7 A plot of the surface V(x, y, z) = h e 



xq = 0.2, i/Q = 0,p X Q = and hi = 1, while in all cases 
the initial value of p y o is always found from the energy 
integral. We see a U type spectrum indicating regular 
motion while the corresponding L.C.E tends asymptot- 
ically to zero. Figure 4a-b is similar to Fig. 3a-b for an 
orbit with initial conditions: xq = 1.0, yo = Q,p x o = 
and h 2 — 2. Here we see a spectrum with a lot of asym- 
metric small and large peaks indicating chaotic motion. 
The plot of the L.C.E vs time in Fig. 4b verifies the re- 
sults of the spectrum. 

The 5(c) spectrum is useful in order to detect is- 
landic and sticky motion. The results are shown in Fig- 
ure 5a-d. Figure 5a shows an orbit with initial condi- 
tions: x — — 1.682, j/o = QiPxO — an d hi — 2. This 
orbit produces the two sets of five islands shown near 
the a;-axis in Fig. lb. Note that the corresponding spec- 
trum produces ten individual regular spectra each one 
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corresponding to an island. As the orbit produces two 
different sets of five islands on the positive and negative 
side of the £-axis, we observe two sets of five spectra 
almost symmetric to the origin. Figure 5b shows the 
spectrum of a nearby orbit. The initial conditions are: 
x a = —1.705, 2/0 = 0,Pa;o = and h 2 =2. This orbit pro- 
duces the two sticky regions around each set of the five 
islands near the x-axis shown in Fig. lb. As there are 
two areas of sticky motion we see two separate sticky 
spectra with asymmetric large and small peaks. Here 
the time is T = 1000 time units. Figure 5c is similar 
to 5b but when T = 4200 time units. We see two dif- 
ferent asymmetric spectra. This indicates that the test 
particle is still in the sticky region. Figure 5d shows 
the spectrum when T = 4300 time units. In this case, 
we see that the two spectra have joined to produce a 
unified chaotic spectrum. This indicates that the test 
particle has left the sticky region and has gone into the 
chaotic sea. The L.C.E vs time of the orbit described 
above is shown in Figure 6. All L.C.Es were computed 
for a time period of 10 5 time units. 

The above results strongly suggest that the S(c) 
spectrum is a powerful dynamical parameter providing 
useful and reliable results for distinguishing between 
order and chaos. Furthermore, the S(c) spectrum de- 
scribes very satisfactorily islandic and sticky motion. 



3 The nature of motion in the 3D model 

In this Section we shall investigate the regular or chaotic 
behavior of the orbits in the 3D potential (1). Working 
in a similar way as in Caranicolas and Varvoglis [9], we 
find that this potential has an energy of escape given 
by 

h esc = 2e + ^ + v/2e (2e + ^) + ^ ^+Z. ( 8 ) 

For values of h > h esc the zero velocity surface is open 
and the test particle is free to escape to infinity. In this 
paper we investigate only bounded motion that is we 
consider always h < h esc . For the adopted values of the 
parameters cj — l,e — 1 we find that h esc — 7.174235. 
Figure 7 shows the surface V(x, y, z) — h esc . Note that 
in the potential (1) the straight lines x = ±y = ±z are 
exact periodic orbits. Using the origin x = y = z = 
as a starting point wc find from Hamiltonian (3) that 
the starting value of the corresponding velocities are 

1/2 

PxQ = Pzo = (x) ' wnue the value of p yQ is obtained 
from the energy integral. 

In order to find the regular or chaotic nature of the 
3D orbits we choose initial conditions (xo,p x o, z o),yo = 



Pzo = 0, where (xo,p x o) is a point on the phase plane of 
the 2D system. This point lies inside the limiting curve 

1 - P l + V{x) = h 2 , (9) 

which is the curve containing all the invariant curves of 
the 2D system. We use h = h 2 and the value of p y o, for 
all orbits, is obtained from the energy integral (3). The 
numerical calculations suggest, that orbits with initial 
conditions {x ,p x0 , z ), y = p z o = 0, such as (x ,p x0 ) 
is a point in the chaotic regions of Fig. 1 a-d and for all 
permissible values of z , produce chaotic orbits. 

On the other hand, it would be also interesting to in- 
vestigate the character of orbits with initial conditions 
{x ,Px0, za),Vo = Pzo = 0, such as (x ,p x o) is a point in 
the regular regions of Fig. 1 a-d. We shall use the S(c) 
spectrum and in order to verify the effectiveness of the 
S(c) spectrum in 3D dynamical systems, we shall com- 
pare the results with two other indicators, the classical 
method of the L.C.E and the spectral method P(f), 
used by Karanis and Vozikis [18]. This method, uses 
the Fast Fourier Transform (FFT) of a series of time 
intervals, each one representing the time that elapsed 
between two successive points on the Poincare (x,p x ) 
phase plane, for 2D systems, while for 3D systems they 
take two successive points on the plane z = 0. Note that 
the coupling of the third component, z, carrying all the 
information about the 3D orbits, is hidden in the def- 
inition of the S(c) spectrum, but in any case it affects 
the values of x,p x and p y entering Eq. (6). Thus, the 
S(c) spectrum provides implicit results for the 3D or- 
bits. On the other hand, using the new S(w) spectrum 
the outcomes are explicit, since the coupling of the third 
component z is located directly in its definition. 

Figure 8 a-d shows (a) the 3D orbit, (b) the cor- 
responding L.C.E. (c) the S(c) spectrum and (d) the 
P(f) indicator, of an orbit with initial conditions: x = 
0.1, p^o = 0, z = 0.1, i/o = p z0 = 0. The value of energy 
is h = 1 , while in all cases the value of p y o is found from 
the energy integral. All dynamical indicators coincide 
that this orbit is regular. Figure 9a-d is similar to Fig. 
8a-d but for a chaotic 3D orbit. Initial conditions are: 
x a = 0.75, p x o = 0, z = 0.1,7/o = Pzo = 0, h = 1. Here 
one observes an asymmetric S(c) spectrum with a lot 
of small and large peaks indicating chaotic motion. The 
L.C.E in Figure 9b and the P(f) indicator in Figure 9d, 
coincide with the results of the spectrum. 

Finally, we shall study the character of orbits start- 
ing near the 3D exact periodic orbits. Figure 10 a- 
d is similar to Fig. 8 a-d but for an orbit starting 
near the 3D exact periodic orbit. Initial conditions are: 

1/2 

%o = 2/o = z = 0,p x0 = (^) ,Pz0 = Pxo + 0.01, 
h = 1. The L.C.E, the S(c) spectrum and the P(f) 
indicator suggest for regular motion. Figure 11 a-d is 
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similar to Fig. 9 a-d but when h = 5. Here the dynam- 
ical indicators reveal the chaotic nature of this orbit. 
The numerical calculations indicate that orbits start- 



ing near the exact periodic orbits become chaotic when 
h > 4. Figure 12 shows a plot of the L.C.E vs h for the 
chaotic orbits in the 3D system. We see that the L.C.E 



8 



Nicolaos D. Caranicolas & Euaggelos E. Zotos 




Fig. 10 (a-d): Similar to Fig. 8a-d, but for an orbit starting near the exact periodic orbit. The motion is regular for h = 1. 




Fig. 11 (a-d): Similar to Fig. 9a-d, but for an orbit starting near the exact periodic orbit. The motion is chaotic for h = 5. 

increases linearly as h increases. Note that the values 4 Discussion and conclusions 
of the L.C.Es of the 3D system are smaller than the 
corresponding values of the L.C.Es of the 2D system. 

Potentials composed of harmonic oscillators have been 
frequently used during the last years by several researchers 
[4, 12-16]. In this article we have studied the regular or 
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chaotic character of orbits in a 3D potential made up of 
perturbed elliptic oscillators. This potential displays ex- 
act periodic orbits, that are straight lines going through 
the origin, together with extended chaotic regions. 

Our investigation begins from the 2D potential. In 
order to keep things simple we tried to study the regu- 
lar or chaotic character of orbits for fixed values of the 
two other parameters and different values of the en- 
ergy. Our numerical calculations have shown that the 
extend of chaotic regions increases rapidly as the value 
of the energy increases. On the other hand, we have 
found that the degree of chaos also increases as the en- 
ergy of increases. The numerical experiments show that 
the L.C.E increases linearly with /i 2 . Furthermore, the 
application of the S(c) spectrum shows that this dy- 
namical parameter is a fast and reliable indicator in 
order to distinguish order from chaos. Note that in or- 
der to obtain the S(c) spectrum we need only one orbit 
and 3000-10000 time units, while in order to obtain the 
L.C.E we use two orbits and much longer time periods. 
Moreover, the 5(c) spectrum describes in a very satis- 
factorily way islandic and sticky motion. On the other 
hand, the time needed in order to obtain the P(f) in- 
dicator, is of order of 5000 time units. 

Using the results of the 2D potential we proceeded 
to study the 3D model. First we tried to find the char- 
acter of the 3D orbits with h = hi and initial conditions 
(ar ,PxO,zo),2/o = Pzo = 0, such as (x ,p x0 ) is a point in 
(i) in the chaotic regions or (ii) in the regular regions of 
Fig. 1 a-d. The numerical results indicate that all orbits 
of type (i) are chaotic. The orbits of type (ii) display 
a different behavior. The numerical outcomes suggest 
that for small values of z (z < 0.1) the orbits remain 
regular, while for larger values of z (z > 0.1) the motion 
becomes chaotic. 

The regular or chaotic nature of the 3D orbits can 
be found using the S(c) spectrum. Comparison with 
the L.C.E and the P(f) spectral method, shows that 
the S(c) spectrum is a good and fast indicator in or- 
der to distinguish regular from chaotic motion, in 3D 
dynamical systems. In the present research, we have 
choose to use dynamical parameters, such as the S(c) 
spectrum and the P(f) indicator, in order to charac- 
terize the nature of orbits, because the main advantage 
of these new methods is that they can provide reliable 
and conclusive results in short integration time peri- 
ods, that is 10 3 — 10 4 time units. Here, we must point 
out that there are also other new dynamical indicators, 
such as the FLI [17], SALI [21] and GALI [22]. These 
indicators need at most 500 to 1000 time units of in- 
tegration time, in order to reveal the true nature of 
an orbit in Hamiltonian systems. The main disadvan- 
tage of the above methods is that we can not use them 



0.5 



0.4 



0.2 



0.1 [. 

1 2 3 4 5 6 7 

h 

Fig. 12 A plot of L.C.E vs h. 

to identify resonant orbits of higher multiplicity which 
correspond to multiple islands of invariant curves in the 
phase plane. Moreover, FLI and GALI are not very sen- 
sitive in the case of sticky motion. On the contrary, the 
SALI method can provide reliable information about 
sticky orbits, using shorter integration time than the 
S(c) spectrum. 

Note that an interesting characteristic of this inves- 
tigations is that both, the 2D and 3D potentials display 
exact periodic orbits and large chaotic regions. The nu- 
merical results, in the 2D and the 3D model, show that, 
for small values of the energy, orbits starting near the 
exact periodic orbits remain regular, while for larger 
values of the energy become chaotic. 
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